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Abstract 

We study the effect of inhomogeneities on light propagation. The 
Sachs equations are solved numerically in the Swiss-Cheese models 
with inhomogeneities modelled by the Lemaitre-Tolman solutions. Our 
results imply that, within the models we study, inhomogeneities may 
partially mimic the accelerated expansion of the Universe provided the 
light propagates through regions with lower than the average density. 
The effect of inhomogeneities is small and full randomization of the 
photons' trajectories reduces it to an insignificant level. 

1 Introduction 

The observations of type la supernovae interpreted within isotropic and 
homogeneous cosmological models imply the accelerated expansion of the 
Universe. These observations \20\ [T7] were made in late nineties and ig- 
nited a "megabit bomb"^ of an enormous number of publications. Among 
many hypothesis that were proposed, the most conservative one assumed 
that inhomogeneities in the energy distribution may mimic the accelerated 
expansion. Although such a hypothesis does not solve the old cosmological 
constant problem (why this constant is so small), it gives a hope for under- 
standing why the vacuum energy density in the concordance model is of the 
same order as the present matter energy density |19j . 

Inhomogeneities may have a twofold effect. Firstly, the averaging pro- 
cedure in general relativity is not well understood yet. Hence, assuming 
homogeneity and then solving the Einstein equations could not lead to the 
proper metric |11]. Secondly, the light propagates differently in inhomo- 
geneous spacetimes. This may modify the luminosity distance - redshift 
relation that is crucial for an interpretation of the type la supernovae data. 



^Stanislaw Lem, Summa technologiae, 1964. 



In this article, we follow [lOlIISlliaEllIlEHlElEl] and study exact 
solutions to the Einstein equations, the so-called Swiss-Cheese (SC) models. 
In such models, inhomogeneities do not influence the global dynamics by 
construction. Therefore, the averaging problem will not be investigated here. 
The SC models provide convenient settings for studies of light propagation 
in inhomogeneous spacetimes. 

The SC models are constructed out of the Einstein-de Sitter (EdS) solu- 
tion with spherical regions removed. The Lemaitre-Tolman (LT) solutions 
are matched to the spacetime in the excised regions.^ Since inhomogeneities 
in the real Universe are not spherically symmetric, it is not obvious how to 
choose density profiles of the inhomogeneous regions. Therefore, we treat 
the SC model as a toy model of the Universe and we search for a reasonable 
"extremal" setting to determine the maximal effect of inhomogeneities on 
the luminosity distance - redshift curve. If shear is neglected, the upper 
bound on the luminosity distance for a given redshift is determined by the 
so-called empty 6eam formula 

The numerical code that we have developed give us large freedom in the 
construction of models. The light may travel non-radially trough arbitrary 
size inhomogeneous regions whose centres do not have to lie in a plane or 
on a regular lattice. This allows us to investigate more general settings 
than these presented so far in the literature. We solve numerically the fully 
relativistic system of equations. 



New exact solutions to the Einstein equations may be constructed out of 
the old ones with a help of a gluing technique. In this article, we consider 
inhomogeneous cosmological models that are made of the LT solutions. The 
large scale evolution is given by the EdS solution that belongs to the LT 
class. Therefore, we start with a short description of the LT solutions and 
matching conditions within this class. 

The LT solutions are spherically symmetric solutions of the Einstein 
equations with a dust source |14^ \23\ Hj. The corresponding line element 
takes the following form in comoving coordinates 



2 Model 



ds^ = dt^ 






1 + 2E 



^For another possibility see 

^For light bundles which have not passed through a caustic [22| . 
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E = E{r) is an arbitrary function and R = R{t, r) satisfies equations 

Rl = 2E+^-f-\KR\ (2) 

''"^ = /Pit' 

where M = M{r) is one more arbitrary function, A is the cosmological 
constant and e = e{t, r) is a dust mass density. The solutions of Q are given 
up to a third arbitrary function t_B(r). Moreover, ([T]), ^ are covariant under 
the coordinate transformation r = f{r'). In order to obtain a particular 
solution one has to specify three functions E(r), M{r) and tBir). Also other 
possibilities exist, e.g. in this article we will set E{r), e{tinit, r) and R{tinit, f) 
at some tinit- All these functions have a simple physical interpretation [Hll8j. 
In principle one may specify them freely, but a general choice will lead to 
pathologies. Supplementary conditions can assure regularity at the centre 
and exclude shell-crossing singularities [18j. 

For A = 0, the solutions of ^ exist in an explicit form |18j . 
When E{r) < 

, , M , , 

R[t,r) = -^(1-cos??), 

(-2^)3/2 

ij-sinrj = ^—^^{t-tB{r)). (4) 



If E{r) = 0, then 



When E{r) > 



1/3 

Rit,r)={"-M{r){t-tBir)r] . (5) 



R{t,r) = ^(coshr?-!), 

(2£;)3/2 

smh7?-7/ = — • 

Two LT solutions can be joined smoothly on the spherical hypersurface 
S given by r = r^. The Darmois junction conditions ^ state that the first 
fundamental forms (intrinsic metrics) and the second fundamental forms 
(extrinsic curvatures) calculated in terms of the coordinates on S should be 
the same on both sides of S. One may show that the matching conditions 
reduce to 

^i|e = ^2|e, Mi|e = M2|s, {tshk = {tB)2k , (7) 
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where indices 1, 2 number matched solutions. The evolution outside the 
matching surface S does not depend on the evolution inside S. The matching 
may be repeated arbitrary number of times as long as different S's do not 
overlap. 

In this article, we study the EdS solution with spherical regions excised. 
Inside each excised region the non-homogeneous LT solution is matched. 
Since the EdS spacetime is homogeneous, the matching surfaces do not have 
to be centred at the same point, but they may be scattered "like holes in 
a cheese". This construction leads to the non-homogeneous model of the 
Universe that is called the Swiss-Cheese (SC) model. 



2.1 Setting 

Hereafter, we assume A = 0. The large scale evolution is given by the LT 
solution with 

E{r) = 0, e{t,r) = ^, t^W = . (8) 

This choice corresponds to the EdS model. At some time tinu we match 
inhomogeneities. They are modelled by the following LT solution. We make 
the same choice as in jl5] and take 

^ 9 

1 I 



e{tinit,r) = Ae^p H ^ 7 +5 . (9) 

Moreover, we assume R(tinii^T^ — r and this together with ( |6j ) determines 
tBir). The formula for E{r) follows from the assumption that at i = tinit 
the speed of the angular expansion of the inhomogeneous regions equals to 
the speed of the expansion of the homogeneous part. It is 

Summarizing, the profile of inhomogeneities at time Unit is given in comoving 
coordinates by the size of the inhomogeneity r^, the position of a peak of 
the density rjv/, it's amplitude A and the width <;", and finally, the parameter 
6 that controls the density at the centres of the inhomogeneous regions. 
The parameters A, rjy.f are not arbitrary, but they are chosen to satisfy the 
matching conditions ([T]). We have found 

^= — ^-5 expH . ) , (11) 
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where rM is determined by solving numerically the equation Mi|s = M2\s- 
This equation may be written with a help of an error function as 



/ 2 



-3e%^^/2^.(ri, + -')iErfC-^) + ^-/(^))) = • (12) 

Mils was found by integration of ([s]) with an assumption M(0) = and 
-^2 Is = §73" is a standard formula for the EdS spacetime. 

init 

We cover the spacetime with many spherically symmetric comoving co- 
ordinate patches of the size centred at inhomogeneities. In principle, the 
free parameters r^, ri,, <j, 6 may vary from one inhomogeneity to another. 
However, in order to reduce the number of degrees of freedom, we assume 
that for a given version of our model the ratio ra/r^ = Sa is fixed. The ratio 
<j/rb = Sc; and S are chosen to be the same for all models studied in this 
article. We do not define the whole spacetime in advance, but construct 
it along the light beam (the congruence of null geodesies) that is evolved 
backward in time (from an observation to an emission). Therefore, for each 
inhomogeneous region there are two additional parameters. One of them, 
plays a role of an impact parameter of the light beam (c^ is a normalized^ 
impact parameter, such that = 1 corresponds to the maximum value of 
the impact parameter). The second parameter is a phase of shear (argcr) of 
the beam at matching surfaces. A change of argcr correspond to a rotation 
of the principal axes of shear in Sachs basis. If argcr is not continuous, 
then centres of inhomogeneous regions may not lie on a two dimensional 
surface. The both parameters c^, argcr are crucial, together with r^, Sa, for 
the relative positions of inhomogeneous regions in spacetime. Finally, it is 
necessary to specify a present moment to (or equivalently Hq = ^ — the 
Hubble constant outside inhomogeneities at the present moment). 

The choice of numerical values of our free parameters is not completely 
arbitrary. We would like to avoid shell crossings in inhomogeneous regions 
and we do not want to have focal points for the congruence of null geodesies 
that we are going to study. In addition, the set of parameters that we choose 
should not be in obvious contradiction to observational data. The allowed 
values of parameters may be determined by the trial and error method. 
However, in order to be able to compare our results to [T5] we make similar 



At each entry to the new coordinate patch = (1 + z)R{t,ra)c^, where wiU be 



defined in Section 3.1 
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choice and the LT inhomogeneous solutions describe central large underden- 
sities surrounded by overdense shells. Of course, there are many other much 
different "good" settings. 

Let us start with choosing units such that the numerical value of the dust 
density at the time of the observation to equals to 1. Thus, we set without 
loss of generality^ Iq = 1/^/6tt. We assume that tinu = 0.2tQ. Next, we set 
s,j = 0.1, 5 = 0.0025 and the matching gives fixed ratio rM/fb = 0.881258. 
In our work, we consider two values of the Sa parameter: Sa = 1 or Sq = 
1.19048. 

Let N denote the number of inhomogeneous regions. There remain 3A^ 
free parameters ((rb)j=i..,Ar, {c^)i=i...N (arg cT)j=i...7v) in the model. In order 
to reduce this number, we fix size of inhomogeneous regions {rh)i for a given 
version of the model and assume that impact parameters (c^)i=i...Ar are dis- 
tributed in the interval [— C0,c^] with the probability distribution function^ 
|c^/c^| or with the uniform probability density.'' Moreover, the phase of 
shear (arg(T)i=i...Ar is equal to zero or uniformly distributed in the interval 
(— TT, tt]. We assumed that in each coordinate patch the beam should enter 
the inhomogeneous region. Hence, 

= min{rf,/ra, 2n/ra\/l - {n/raY} , (13) 

where the second value prevents inhomogeneities from overlapping and /r^ = 
1/sa- In our setting, we have c*^ = r^/ra ^ 0.84. 

In short, our models were defined in this Section. The differences between 
particular settings studied in this article are presented in Appendix |A| 

^In the EdS model p = 
The so-called "non-aligned" version of our SC models. Inhomogeneities are randomly 
spread in the spacetime and the probability that the beam will enter inhomogeneous region 
with a particular impact parameter is proportional to the circumference 2tic^. Such 
the probability distribution functions provides proper randomization of the path of the 
beam. 

'^The so-called "aligned" version of our SC models. Inhomogeneities are randomly 
spread in the plane in which the beam propagates. The beam is randomized only in the 
plane. 
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3 Luminosity distance and angular diameter dis- 
tance 



The corrected luminosity distance [22j from the source that is moving with 
the four-velocity u'g to the observer with u'q is defined by 

where SAq is the area of a cross-section of the beam seen by the observer (it 
does not depend on Uq) and 5Vts is the solid angle of the beam that is mea- 
sured at the source in a tangent 3-space orthogonal to {5^s depends on 
tig ). Therefore, di depends only on the position of the source, the position 
of the observer and the four-velocity of the source. Similarly, interchanging 
roles of the source and the observer one can define the angular diameter 
distance dA = {^^s / ^^oY^"^ ■ This quantity is not the observer indepen- 
dent {5Vto depends otiuq). The reciprocity theorem [22] proves that in any 
spacetime di = {1 + z)dA, where z is a redshift of the source seen by the 
observer. The uncorrected luminosity distance is given by di = {1 + z)dL- 
The additional redshift factor corresponds to the change of the energy of 
photons. Thus, d^ has physical meaning and may be used to calculate the 
apparent brightness. Now, we may write 



<<.^(l..)(^)\ (.5) 

or using the reciprocity theorem 

In this article, we assume that the observer see the source within some small 
solid angle SVIq and we trace back the beam to the source, evolve A along the 



path and calculate dA- Finally, one may use (16) and the geodesic equation 
to determine di as a function of the redshift z. However, we have decided 
to present results in terms of dA- 

The evolution of the area of the cross-section of the beam is given by the 
Sachs optical equation |21] 



+ VI = 0, (17) 
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where A is an affine parameter, is a wave vector of the beam, a is complex 
shear 

■ \{Vah){V^kP)-9\ (18) 



and 6 is an expansion 



1 



(19) 



We are interested in the evolution backward in time, so the initial conditions 
we adopt are 



'A\o 
dVA 



0, 



(20) 



dX 



o 



o 



where the last equation follows from the definition of a solid angle. The 
observation corresponds to A = 0. The angular diameter distance does not 
depend on 5flo as long as 6Qo is small enough. 



One procedure for solving (17) would be to define initial conditions for 



a congruence of null geodesies (in agreement with (20)), solve the geodesic 



equation (28) with these initial conditions, calculate the expansion (19), 



shear (18) and finally solve (17) with initial conditions given by (20). How- 



ever, it is more convenient to adopt a different approach. Namely, one solves 



the geodesic equation (28) for a single central null geodesic and calculates 



shear and the expansion of the congruence from the remaining Sachs optical 
equations^ 



1 



2 RfjLv^^k 



„2 I i2 

— + r+cjr = 

dX 

^ + 2ea = C^p^sL'^k^WLi 
(iX 



(21) 
(22) 



where Li is the spacelike vector orthogonal to the light ray (one of the 
vectors of the Sachs basis), |o"|a=o = and we set the phase of a to at 
A = 0. In our case, (22) took a particularly simple form because of spherical 



symmetry of the LT solutions. It can be verified that 



(23) 



*It is sufficient to use the equation 6 
assumption (S|a|)|A=o = 0. 



dA/d\/{2A) with (|22| supplemented by the 
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is a good choice of Li (see |5J). In the parametrization ([T]) the right hand 
sides of (21), (22) take the form 



^Rf.ukf'k^ = -iTT{l + zfe, (24) 



2 



C^f.^sL'^k^k'TLi = 2.A'(e-^), (25) 



R ' V 47ri?3 J 

where c,^ corresponds to the value of the impact parameter in a given coor- 
dinate patch. 

The Sachs equations depend on the wave vector of the central light ray in 
the beam, so it is necessary to find the ray trajectory. In the next subsection, 
we describe the geodesic equation. 

In a general spacetime, solving the Sachs and the geodesies equations 
involves numerical calculations. For spacetimes that are "on-average-RW" 
models, the effective procedure was suggested [26l [9]. This procedure ne- 
glects shear and for an "on-average-EdS" model gives the formula 



1 [i+zf-ji+zr^ 

dA{z) = 77 75 > (26) 

^0 2/3(1 + z)4 

where /3 = |\/25 — 24a and a is a dimensionless smoothness parameter (a 
mass- fraction of the matter in the Universe that is not bounded in galaxies) . 
This is the so-called partially filled beam approximation and for a = it is 



known as the empty beam approximation. We will compare (26) to our 
numerical results in Subection 14. 5[ 



3.1 Geodesic equation 

In order to determine null geodesies we will make use of Killing symmetries of 
the spacetime ([l| . It follows from the form of the line element ([!]) that one of 
Killing vectors implied by spherical symmetry takes the form rj^ = {d/d(j))^. 
Let be the four-velocity of a photon and let A be an affine parameter 
along photon trajectory 

_ ( d\ dr fd\ d^ f_d\ dc^ f d_\ 

""'dxydtj^dxydrj^dxydej^dxyd^) ' ^ 

It satisfies the geodesic equation u^V^u'^ = 0. The quantity u^r?^ = 
dcj) / dXR^ siv? 9 = C0 is conserved along the photon path. We can easily 

^This is the effective formula and no rigorous derivation is known at present. 
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choose our coordinate system in such a way that at some point 0{Xp) = 0. 
This imphes that = along the whole trajectory and ^(A) = or 
0(A) = const. Therefore, we assume without loss of generality that one 
of the angular variables is constant along the photon path. We choose our 
coordinate system to have 9 = 7r/2. Thus, d(j)/dX = c^/R^ . Moreover, is 
a null vector and this gives us the first integral of motion u'^u^ = 0. Using 
these equations and t-component of the geodesic equation we obtain the 
system of the first order differential equations 



dt 
dX 



z + l, (28) 



dr _ I 
dX ~ R 



[1 + 2E) (z + l) 



1 
R^ 



de 

dX 







d^ _ C0 

dx ~ W 

2 

dz _ R,rt \^ , ^'f> I -^^rt R,t 

I i + 2 ) + 



dX R^r R^ V R,i' R 

where ^(A) is a new auxiliary function chosen up to an additive constant. 
The second equation is inconvenient for numerical calculations because the 
expression under the square root can be a source of numerical problems for 
near zero values. Moreover, the plus and the minus sign can be encoded in 
initial conditions. Therefore, we use the second order equation for r(A). It 
has the form 

d'r Rrtil + zldr 1 ( .^^.R,rr j,\(,,..2 
dA^ + ^ R, dA + if V^^ + ^^^iJ^-^'^j 1^^" + ^) 

-(1 + 2^)^=0. (29) 

Let be the four- velocity of the source that is comoving with a dust. 
Since there are no off-diagonal terms in the metric ([T]), we have 1 + z{X) = 
u^u^ = w(A), where w(A) is the frequency of a light signal of wave vector u^^ 
measured by the observer comoving with the source and a dust. If z(0) = 0, 
then z{X) can be interpreted as the redshift of the signal that was emitted 
at some A < (by the source comoving with a dust) and was measured by 
the dust comoving observer at A = 0. 
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In this article, we assumed that the observer is at the matching surface 
between the EdS solution and the inhomogeneous region, hence we set r(0) = 
Tq. We have dr/dA < at A = and the remaining initial conditions are: 

m = to, e{o) = f, 0(0) = 0, z{o) = o. 
4 Numerical results 

As a starting point we take the model with five inhomogeneous regions that 
was investigated by Marra, Kolb, Matarrese, Riotto (MKMR) in [15] . It 
was shown there that the MKMR model almost reproduce the angular di- 
ameter distance - redshift curve of the Robertson- Walker (RW) model with 
Qm = 0.6 and Q\ = 0.4. This results suggest that the effect of inhomo- 
geneities may be significant. We would like to verify how it depends on the 
assumptions made and the particular setting used by MKMR. Therefore, 
we generalize the MKMR model step by step and follow the changes in the 
angular diameter distance - redshift curve. The similar analysis was done in 
[25j , where the large effect in the MKMR model was explained as a result of 
insufficient randomization. The main difference between our study and [25| 
comes from the fact that in [25] the weak field gravitational lensing theory 
was used and shear was neglected. We solve numerically the fully relativistic 
system of equations and evaluate directly the effect of shear. 

The differences between particular settings studied in this Section were 
summarized in Appendix [Xj The order of the keys in the figures corresponds 
to the order of the curves for the redshift z = 1.8. 

4.1 MKMR model 

In our setting the MKMR model corresponds to = 0.1, Sa = 1, r^, = 0.042 
and arg cr = = in each inhomogeneous region. We reproduced the 
result of [15J. The change in the angular diameter distance, AdA{z) = 
dA{z) — dA{z)\Rw, compared to the RW model with = 0.6 and Q\ = 0.4 
is presented in Fig. [l} It coincides with Fig. 13 in [15]. 

4.2 Non-radial beams 

In the MKMR model, the light propagates radially in each coordinate patch. 
Since the mass is concentrated in spherical shells and the density in the cen- 
tral part of the inhomogeneity is tiny, the beam travels through regions with 
lower density than the average. The angular diameter distance in such model 
should deviate from the EdS value as it already follows from the weak lens- 
ing analysis. It seems interesting to calculate the angular diameter distance 
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for more typical beams, i.e. the beams that are not so much statisticaUy 
distinguished. Such analysis requires the study of non-radial beams and 
randomization of the sequence of impact parameters. This extension of the 
MKMR is not a trivial one and leads to some subtle problems that will be 
addressed in this Subsection. 

Firstly, for non-radial beams an impact parameter takes random values 
in each inhomogeneous region and centres of those regions are not lined up 
any more. Clearly, gaps between inhomogeneities and an upper limit for 
an impact parameter are necessary to avoid overlapping of inhomogeneities. 
Gaps are introduced by setting Sa = 1.19048 which implies = 0.05. In 
addition, we assume that the absolute value of the impact parameter is not 



larger than = 0.84 (as defined in the equation (13)). This assumption 
improves statistical properties of our models because it assures that the 
beam enters an inhomogeneous region in each coordinate patch. 

Secondly, one has to decide how inhomogeneities are distributed in the 
spacetime. This is defined by the statistical properties of the sequence of 
impact parameters. We consider two possibilities in our paper. The spatial 
average density of the t = const hypersurfaces does not deviate much from 
the EdS value. Hence, the natural assumption is that the average density 
along the beam which goes through inhomogeneous regions should not de- 
viate from the average density along the beam that propagates in the EdS 
spacetime. To achieve this, impact parameters are distributed in [— c^,c^] 
with the probability distribution function |c^/c^|. Such the probability dis- 
tribution function corresponds to the trajectory of the typical beam that 
travels through inhomogeneous regions which are spread randomly in the 
spacetime. The second possibility we consider is that centres of all in- 
homogeneous regions are spread randomly in a plane and that the beam 
propagates in this plane. We will call this version of our models "aligned" 
and we will refer to remaining our models as to "non-aligned" . It is assumed 
that the typical beam in the plane with aligned inhomogeneities corresponds 
to the sequence of impact parameters uniformly distributed in the interval 
[— c^,c^]. In our setting, this implies that the light propagates through re- 
gions with lower average density than the spatial average density in the 
model. Such configuration of inhomogeneities is unnatural, but it is conve- 
nient to model the selection effect, i.e. the light from the sources that are 



^"^We acknowledge private communication with Syksy Rasanen and Krzysztof Bolejko. 
See also [5]. 

^^The finite size of inhomogeneities and the assumption that the beam should enter an 
inhomogeneous region in each coordinate patch imply that the distribution of centres of 
inhomogeneities in spacetime is not uniform. 
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Figure 1: The difference of the angular diameter distances (times Hq) 
between the MKMR model, the radial MKMR model with gaps and the 
RW model with = 0.6, Qa = 0.4, (AdAiz) = dA{z) - 

observed is more likely to travel through low density regions. 

Next, the problem that is sometimes not handled properly in the lit- 
erature is the value of shear at the boundary between coordinate patches. 
In general, the principal axes of shear do not have to coincide with our 
choice of Sachs basis. The rotation of principal axes of shear is necessary in 
non-aligned models and may be realised by a random change of the phase of 
shear (argu) with the uniform probability distribution in (— tt, vr]. In aligned 
models, one should assume that the phase of shear do not change between 
coordinate patches. For a sake of curiosity, we will also consider the third 
possibility: an aligned model (a vanishing phase of shear) with the impact 
parameters probability distribution function |c0/c^|. 

In the remaining part of the paper, the following terminology will be 
used. Whenever we will refer to the aligned version of the non-radial SC 
model, we will explicitly indicate that. We will also explicitly indicate when- 
ever only radial beams will be studied in the model (if different than the 
MKMR model). The remaining models are assumed to be non-aligned and 
to contain non-radial beams. 

Let us start with applying the extension presented in this Subsection to 
the MKMR model. In the first step we add gaps. The angular diameter 
- redshift curve in the radial MKMR model with gaps does not reveal big 
changes — Fig. [TJ The inhomogeneities appear for a little bigger redshift, as 
expected. Next, we randomize impact parameters. The non-radial beams in 
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Figure 2: The typical non-radial trajectories (comoving coordinates) in the 
MKMR model with gaps. The inner circle corresponds to r = rf, and the 
outer to r = Tq. All inhomogeneous regions were overimposed on the single 
plot and the planes of the trajectories were rotated. 



the MKMR model with gaps are curved. We present them in Fig. [2] where c,^ 
is distributed in [—£^,5^] with the probability density \c^/c*^\ (non-aligned 
setting). The angular diameter distance is lower than in the EdS model 
because the beam spends "more time" in the region of higher density — Fig. 
[Sj The effect of inhomogeneities on the angular diameter distance is opposite 
than for radial beams. However, this is not a typical property of non- 
radial beams, but an effect of week randomization (only five inhomogeneous 
regions) . 

4.3 Small inhomogeneities 

In the MKMR model with gaps (non-radial), the angular diameter distance 
dA{z) depends strongly on the sequence of impact parameters. One should 
average dA{z) over many runs to obtain a reliable result (like in [25j) or 
reduce the size of inhomogeneities to obtain the better statistic in a single 
run. We decided to reduce the size of inhomogeneities 100 times (we set 
= 0. 00042). In the model with small inhomogeneous regions (hereafter, 

^^The diameter of the inhomogeneous region is around 10 Mpc at the moment of an 
observation. 
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Figure 3: The difference of the angular diameter distances (times Hq) 
between the empty beam approximation, the SI models, the MKMR models 
with gaps (radial and non-radial), the RW model with Q,m = 0.3, Q\ = 0.7 
and the EdS model {AdA{z) = ^^(z) - dA{z)\Eds)- HqAcIa for the RW 
model was multiplied by a factor 10~^ to make the comparison more explicit. 
The order of the keys corresponds to the order of the curves for the redshift 
z = 1.8. 



the SI model), the density along the beam is a fast varying function — Fig. 
|4j For radial beams, the angular diameter distance does not change much, 
but in a general case the effect of inhomogeneities on the angular diameter 
distance is negligible — Figs [3j [5] For non-radial beams in the aligned SI 
model, the effect of inhomogeneities on light propagation for z < 1.5 is lower 
than 10% of what is needed to explain the accelerated expansion without 
introducing the cosmological constant. It is around one-third of the effect 
in the original MKMR model. Since this is the maximal effect we have 
observed in our models (after randomization of the trajectory of the beam), 
the aligned SI model corresponds to our "extremal" setting. The empty 
beam formula (26) gives the largest angular diameter distance. 

The typical trajectories in the SI model (non-aligned and aligned) are 
presented in Fig. |6j It follows from Figs |3j [5] that the effect of inhomo- 
geneities on the angular diameter distance is reduced in our models to the 
insignificant level by proper randomization of the beam's trajectory. 
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Figure 4: The density along the beam in the model with small inhomoge- 
neous regions (the SI model) and in the model with five large inhomogeneous 
regions (the MKMR model). 
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Figure 5: The difference of the angular diameter distances (times Hq) 
between the SI model and the EdS model, (AdAiz) = dA{z) — dA{z)\Eds)- 
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Figure 6: The typical non-radial trajectories (comoving coordinates) in a) 
the SI model, b) the aligned SI model. The inner circle corresponds to 
r = rfj and the outer to r = r^. The planes of trajectories were rotated and 
all inhomogeneous regions were overimposed on the single plot. 



4.4 Shear 

In the previous subsections, we have introduced non-radial beams with non- 
vanishing shear along the trajectories. In most models studied in the lit- 
erature (including [25j), shear is assumed to be negligible for the angular 
diameter distance - redshift relation. We have verified that this is indeed 
true in our models. 
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Figure 7: The angle between the Sachs basis and the principal axes of shear 
in the SI model. 
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Figure 8: The effect of averaging of shear by inhomogeneities whose centres 
are not aligned. The bottom curve corresponds to the SI model. The middle 
curve corresponds to the aligned SI model. The top curve corresponds to 
the aligned model with the probability distribution function of the impact 
parameters \C(f,/5^\ (the aligned SI* model). 
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Figure 9: The difference of the angular diameter distances (times Hq) 
between the SI model with neglected shear (aligned, non-aligned), the EdS 
model, the aligned and the non-aligned SI model, and the aligned SI* model. 
The order of the keys corresponds to the order of the curves for large red- 
shifts. 

In the non-aligned models, the angle between Sachs basis and principal 
axes of shear is uncorrelated between inhomogeneous regions (Fig. [7]). In 
contrast to that, in the aligned models the phase of shear is constant and 
equal to zero. 

We compare the value of shear in three versions of the SI models: non- 
aligned, aligned, and aligned (a vanishing phase of shear) with the probabil- 
ity distribution function of the impact parameters Icj^/c^j. The last model 
was introduced to verify validity of some calculations presented in the liter- 
ature. It will be distinguished with a symbol SI*. The modulus of shear is 
a smooth function and is presented in Fig. [8| It is interesting to observe in 
this figure the averaging effect of non-aligned inhomogeneities on shear. 

Let n be a number of inhomogeneous regions encountered by the beam. 
At first sight the non-aligned positions of inhomogeneities (in the non- 
aligned SI model) lead to a random walk in |cj| ~ ^/n in contrast to the 
aligned case where |(t| ~ n (the aligned SI model). We have verified that 
the curves in Fig. [8] do not satisfy this relations. Since structures are grow- 
ing in time, the expected value of the change of shear between neighbouring 
inhomogeneities is time dependent. Therefore, this process may be studied 
as a continuous time random walk. 

In Fig. [9] we present the difference of the angular diameter distances 
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Figure 10: The difference of the angular diameter distances (times Hq) in 
four models with different sizes of inhomogeneous regions SI, SI2, 813, SI4 
and the EdS model (Adyi(z) = dA{z) — dAiz)Eds)- The order of the keys 
corresponds to the order of the curves for the redshift z = 1.8. 



between the SI model, the aligned SI model, the aligned SI* model and the 
SI model with neglected shear (non-aligned and aligned), and the EdS model. 
This comparison reveals that the effect of shear for non-aligned SI models 
is tiny (for z = 1.5 it less than 0.1% of the correction to the EdS angular 
diameter distance). For aligned SI models, it is less than 1% (for z = 1.5). 
However, the aligned SI* model (with the probability distribution function 
\C(j,/c*^\), overestimate it around two orders of magnitude (for z = 1.5 the 
effect of shear is around 10% of the correction to the EdS angular diameter 
distance). 

The random number generator in our code was initialized with the same 
random seed. Therefore, three bottom curves in Fig. [9] are not jagged for 
small redshifts. 



4.5 Partially filled beam approximation and diff^erent sizes 
of inhomogeneities 

It was argued in [25] that the dimming of supernovae in the MKMR model 
may be roughly estimated with a help of the partially filled beam approx- 



imation (26). We apply this approximation to five models: the SI model, 



three Sin models with inhomogeneities n-times bigger than in the SI model 
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Figure 11: The partially filled beam formula compared to the numericahy 
calculated angular diameter distance for the SI, SI2, SI3, SI4, MKMR 

models, {Ad'^{z) = d a{z) partially filled beam - dA{z), Ma{z) = dA{z)- 

{dA{z))Eds)- The smoothness parameter was calculated for each model us- 
ing ( 30 ) . The order of the keys corresponds to the order of the curves for 
the redshift z = 1.8 and z = 1.7. 



(where n = 2, 3, 4), and the MKMR model. The angular diameter distance, 
in variations of the SI model, seems to depend slightly on the size of in- 



homogeneities (Fig. 10 ) and the sequence of impact parameters induces big 
statistical fluctuations.^^ We suppose that it is not size of inhomogeneities 
that matters, but randomization, i.e. how the average density along the 
beam deviates from the EdS value. This may be seen as follows. Let us 
calculate the smoothness parameter a for each model. The smoothness pa- 
rameter is a mass- fraction of the matter in the Universe that is not bounded 
in galaxies. We calculate it for our SC models using 

" = ^ / n — mT^^ ' 



^^The code was initialized with the same random seed, but sizes of inhomogeneities are 
different. Therefore, the different parts of a random sequence corresponds to different 



redshifts 
i*The 

as may be seen in Fig. [3] 



'^'^The radial versions of MKMR and Sf models give similar angular diameter distance 



21 



where Aq = denotes an affine parameter at the observer and Ag corre- 
sponds to the emission. In our models, the spatial average density of the 
t = const hypersurfaces is approximately equal to the density in the EdS 
model. Therefore, one may expect that a should not diverge far from 1. 
However, the bigger are inhomogeneities, the weaker is randomization of 
the path of the beam. Since low density regions occupy larger volume in 
the studied models, it is more likely that a will decrease with increasing 
size of inhomogeneities (the variance will grow). We have found a to be 
approximately equal to 0.964, 0.939, 0.915, 0.895, 0.815, 0.412 for the SI, 
SI2, SI3, SI4, ahgned SI, MKMR models, respectively. For SI models that 
contain hundreds of inhomogeneities, a is slightly lower than 1. The values 
of a for particular models together with Fig. 10 suggest that the effect of 
inhomogeneities on the angular diameter distance depends on the average 
density along the beam (as expected from the week lensing approximation). 

We plotted Fig. |11| which shows the difference between the partially filled 
beam formula and the numerical result for the approximated model. It 
reveals that the partially filled beam approximation is acceptable for smaller 
values of a. The accuracy is around 1% of the correction to the EdS angular 
diameter distance for the MKMR model and 8% for the aligned SI model 
(for z = 1.5). In our models, the partially filled beam approximation does 
not work for larger values^^ of a. 



5 Summary 

In this article, the effect of inhomogeneities on light propagation was inves- 
tigated in the framework of the SC models. We have examined the angular 
diameter distance - redshift relation. This type of distance is related to the 
luminosity distance by the reciprocity theorem. Therefore, the theoretical 
angular diameter distance - redshift relation is crucial for an interpretation 
of the type la supernovae data. 

Our analysis confirms that inhomogeneities may partially mimic the ac- 
celerated expansion of the Universe provided the light propagates through 
regions with lower than the average density. The effect is small and it be- 
comes negligible if the average density along the beam does not deviate from 
the corresponding EdS value. In light of our work and the weak field grav- 
itational lensing analysis [25], the result [15] that suggest more significant 
influence of inhomogeneities is due to a peculiar setting of the underlying 

^^For df ~ 1, the effect of inhomogeneities seems to be smaller than statistical fluctua- 
tions. 
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model. The precise size of the effect depends on the details of the model 
that was studied. Since the SC models are toy-models of the real Universe, 
it is speculative to base on them the final conclusions. Nevertheless, what 
our analysis shows is that, within the models we have studied, the effect 
of inhomogeneities remains too small^^ to explain the type la supernovae 
observations without dark energy. Our analysis of the fully relativistic and 
non-linear models did not reveal any stronger effect on the angular diame- 
ter distance than that predicted by the partially filled beam approximation. 
Randomization reduces the effect considerably in accordance with |25) . [7], 
[6] (and the others). 

We have directly evaluated the effect of shear on the angular diameter 
distance. Within our models, the effect of shear was negligible, but the 
models that do not take into account the rotation of the principal axes of 
shear (e.g. |24j ) may overestimate its effect around two orders of magnitude. 
In these models, the overestimated shear plays a minor role and it may lead 
to the small underestimation of the effect of inhomogeneities. 

Our results suggest that the size of inhomogeneities is not crucial for the 
angular diameter distance, provided that non-radial models are sufficiently 
randomized. We have found that the partially filled beam formula ( 26 ) gives 
good approximation to the angular diameter distance if the average density 
along the beam is much lower than the corresponding EdS value. 

Finally, we stress the analysis presented in this article does not touch 
directly the "averaging problem" in general relativity. The SC models be- 
have on large scales by construction as the RW models, thus the influence 
of inhomogeneities on the global expansion rate of the Universe cannot be 
studied within this framework. 
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A Models and parameters 
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Impact parameters were uniformly distributed in [^c^,c^], where — 0.84. 



Table 1: The settings of our models. The remaining parameters coincide 
for ah models and were defined in Section[2| The parameters c^, arg((j) are 
randomly chosen at each entry to a new coordinate patch. If not indicated, 
the default probability distribution of the impact parameter Ctf, is \C(j,/c*^\ 
(in the interval [— c^,c^]), where = 0.84. The phase of shear arg(cj) is 
uniformly distributed in (— 7r,7r]. 
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